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Solutions obtained by the quasilinearization method (QLM) are compared with the WKB solu- 
tions. While the WKB method generates an expansion in powers of h, the quasilinearization method 
(QLM) approaches the solution of the nonlinear equation obtained by casting the Schrodinger equa- 
tion into the Riccati form by approximating nonlinear terms by a sequence of linear ones. It does 
not rely on the existence of any kind of smallness parameter. It also, unlike the WKB, displays no 
unphysical turning point singularities. It is shown that both energies and wave functions obtained 
in the first QLM iteration are accurate to a few parts of the percent. Since the first QLM iterate 
is represented by the closed expression it allows to estimate analytically and precisely the role of 
different parameters, and influence of their variation on the properties of the quantum systems. The 
next iterates display very fast quadratic convergence so that accuracy of energies and wave functions 
obtained after a few iterations is extremely high, reaching 20 significant figures for the energy of the 
sixth iterate. It is therefore demonstrated that the QLM method could be preferable over the usual 
WKB method. 



PACS numbers: 03.65.Ca, 03.65.Go, 03.65.Sq 



I. INTRODUCTION 

The quasilinearization method (QLM) was constructed 
as a generalization of the Newton-Raphson method 0, Q 
for the nonlinear differential equations to yield rapid 
quadratic and often monotonic convergence to the ex- 
act solution. It was developed originally in theory of 
linear programming by Bellman and Kalaba H 0] to 
solve nonlinear ordinary and partial differential equations 
and their systems. In the original works of Bellman and 
Kalaba 0,3, however, the convergence of the method 
has been proven only under rather restrictive conditions 
of small intervals and bounded, nonsingular forces [Tol | 
which generally are not fulfilled in physical applications. 
This could explain an extremely sparse use of the tech- 
nique in physics, where only a few examples of the ref- 
erences to it could be found 

011 011 SI Recently, 
however, it was shown |ldj | by one of the present au- 
thors (VBM) that a different proof of the convergence 
can be provided which allows to extend the applicability 
of the method to realistic forces defined on infinite in- 
tervals with possible singularities at certain points. This 
proof was gen eralized and elaborated in the subsequent 
works [3 El III Ell- 
in the first paper of the series E3j the analytic re- 
sults of the quasilinearization approach were applied to 
the nonlinear Calogcro equation [5j in the variable phase 
approach to quantum mechanics, and the results were 
compared with those of the perturbation theory and with 
the exact solutions. It was shown that the number of the 
exactly reproduced perturbation terms doubles with each 
subsequent QLM iteration, which, of course, is a direct 
consequence of a quadratic convergence. 

The numerical calculation of higher QLM approxima- 
tions to solutions of the Calogero equation with different 
singular and nonsingular, attractive and repulsive poten- 



tials performed in the work |l2| has shown that already 
the first few iterations provide accurate and numerically 
stable answers for any values of the coupling constant 
and that the number of iterations necessary to reach a 
given precision increases only slowly with the coupling 
strength. It was verified that the method provides accu- 
rate and stable answers even for super singular potentials 
for which each term of the perturbation theory diverges 
and the perturbation expansion consequently does not 
exist. 

In the third paper of the series |l-3| the quasilineariza- 
tion method was applied to other well known typical 
nonlinear ordinary differential equations in physics, such 
as the Blasius, Duffing, Lane-Emden and Thomas-Fermi 
equations which have been and still are extensively stud- 
ied in the literature. These equations, unlike the nonlin- 
ear Calogero equation [B| considered in references [Tol[l^| . 
contain not only quadratic nonlinear terms but various 
other forms of nonlinearity and not only the first, but 
also higher derivatives. It was shown that again just a 
small number of the QLM iterations yield fast convergent 
and uniformly excellent and stable numerical results. 

In the work [l4| the quasilinearization method was ap- 
plied to quantum mechanics by casting the Schrodinger 
equation in the nonlinear Riccati form and calculating 
the QLM approximations to bound state energies and 
wave functions for a variety of potentials, most of which 
are not treatable with the help of the perturbation the- 
ory or the 1/N expansion scheme. It was shown that the 
convergence of the QLM expansion for both energies and 
wave functions is very fast and that already the first few 
iterations yield extremely precise results. In addition it 
was verified that the higher QLM approximations, unlike 
those in 1/N expansion method, arc not divergent at any 
order. 

The present work is devoted to comparison of QLM 
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and WKB. Indeed, the derivation of the WKB solution 
starts by casting the radial Schrodinger equation into 
nonlinear Riccati form and solving that equation by ex- 
pansion in powers of h. It is interesting instead to solve 
this nonlinear equation with the help of the quasilin- 
earization technique and compare with the WKB results. 
Such a procedure was performed in works 0, @ , where it 
was shown that the first QLM iteration reproduces the 
structure of the WKB scries generating an infinite series 
of the WKB terms, but with different coefficients. Be- 
sides being a better approximation, the first QLM itera- 
tion is also expressible in a closed integral form. Similar 
conclusions are reached for higher QLM approximations 
and it can be shown ^| that the p-th QLM iteration 
yields the correct structure of the infinite WKB series 
and reproduces 2 P terms of the expansion of the solution 
in powers of h exactly, as well as a similar number of 
terms approximately. 

That the first QLM iteration already provides a much 
better approximation to the exact solution than the usual 
WKB is obvious, not only from comparison of terms of 
the QLM and WKB series HUE!], but also from the fact 
that the quantization condition in the first QLM iteration 
leads to exact energies for many potentials 0, 0] such 
as for the Coulomb, harmonic oscillator, Poschl- Teller, 
Hulthcn, Hyllcraas, Morse, Eckart and some other well 
known physical potentials, which have a simple analytic 
structure. By comparison, the WKB approximation re- 
produces exact energies only in the case of the first two 
potentials. 

The goal of this work is to point out that also for 
other potentials with more complicated analytical struc- 
ture QLM iterates provide much better approximation 
than the usual WKB. If the initial QLM guess is properly 
chosen the wave function in all QLM iterations, unlike the 
WKB wave function, is free of unphysical turning point 
singularities. Since the first QLM iteration is given by 
an analytic expression 0,H GjJ, GH E 13 , it allows one 
to analytically estimate the role of different parameters 
and the influence of their variation on different charac- 
teristics of a quantum system. The next iterates display 
very fast quadratic convergence so that accuracy of en- 
ergies obtained after a few iterations is extremely high, 
reaching up to 20 significant figures for a sixth iterate as 
we show on the example of different widely used physical 
potentials. 

The paper is arranged as follows: in the second chapter 
we present the main features of the quasilinearization ap- 
proach to the solution of the Schrodinger equation, while 
in the third chapter we consider the application of the 
method to computations with the anharmonic oscillator, 
logarithmic, two-power (double- well) , and Wood-Saxon 
potentials and to the two-body Dirac equation with static 
Coulomb potential. The final, forth chapter is devoted 
to the discussion of the results, convergence patterns, nu- 
merical stability, advantages of the method and its pos- 
sible future applications. 



II. QUASILINEARIZATION METHOD 

The usual WKB substitution 

x(r) = Cexp I A / y(r')dr' 



(1) 



converts the Schrodinger equation to nonlinear Riccati 
form 



dy(z) 
dz 



+ (k 2 (z)+y 2 (z))=0. 



(2) 



Here k 2 (z) = E-V -1(1 + 1) /z 2 , X 2 = 2m/h 2 and z = Xr. 

The proper bound state boundary condition for poten- 
tials falling off at z ~ zq ~ oo is y(z) = const at z > zq. 
This means that y'(zo) = 0, so that Eq. J2J at z ~ zq 
reduces to k(zo) 2 + y 2 (z^)) = or y{z )) = ±ik(z ). We 
choose here to define the boundary condition with the 
plus sign, so that y(zo) = i k(zn ). 

The quasilinearization 0,13 El of this equation gives 
a set of recurrence differential equations 



dy P {z 
dz 



y 2 p _ 1 (z)-2y p (z)y p _ 1 {z)-k 2 (z) (3) 



with the boundary condition y p (zo) = ik(zo)- 

The analytic solution of these equations expresses 
the p-th iterate y p (z) in terms of the previous iterate: 



Vp(z) 
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fp-i(z) - 
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(4) 



Indeed, differentiation of both parts of Eq. (0} leads im- 
mediately to Eq. J2J which proves that y p (z) is a solution 
of this equation. The boundary condition is obviously 
satisfied automatically. 

To utilize the recurrence relation for wave function 
computation one has to pick up a proper initial guess. For 
the zeroth iterate yo(z) it seems natural to choose the 
zero WKB approximation that is to set yo(z) = ik(z), 
which in addition automatically satisfies the boundary 
condition. However, one has to be aware that this choice 
has unphysical turning point singularities. According 
to the existence theorem for linear differential equations 
[3, if y p -i(z) in Eq. © is a discontinuous function of 
z in a certain interval, then y p (z) or its derivatives may 
also be discontinuous functions in this interval, so conse- 
quently the turning point singularities of yo(z) may prop- 
agate to the next iterates. To avoid this we choose [l(| 
the Langer WKB wave function as the zero iteration. 
This function near the turning p_oints a and b is given by 
the simple analytic expression 



Xi{r) 



Si(r) = -X 



\k(r)\ 

\k(s)\ds 



dS?(r) 



(5) 
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Here Ai denotes the Airy function, i = a,b, k 2 {r) = 
2m{E - V(r)) - {I + l/2) 2 /r 2 , d is -1 for a < r < b, 
and 1 for r < a, r > 6, and c a = 1, c& = (—1)™, where 
n = 0, 1, 2, . . . is the number of the bound state. Xa{r) 
and Xb(r) are continuous across the turning points and 
coincide with the usual WKB solution far from them. It 
is easy to check that Xair) and Xb( r ) coincide at some 
point in the interval (a, b) between the turning points, 
and that their values, but not derivatives, can be matched 
at that point. 



III. EXAMPLES 

To show that the first QLM iteration gives very ac- 
curate results for wave functions and energies, as well as 
demonstrate very fast convergence of the next iterates let 
us consider five typical examples of potentials of rather 
different form used in atomic, nuclear and quark physics. 

Let us start from the anharmonic oscillator V(r) = 

1„5 




FIG. 1: Comparison of the Langer WKB solution xo (dashed 
curve), the exact solution x exact (dotted curve) and the first 
QLM iterate Xm„ (solid curve) for the ground state of the 
anharmonic oscillator V(r) = \r 5 . The last two are indistin- 
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This potential is typically used in different nu- gu i s hable on the plot. Here x = nr, k 2 = 2mE/h 
clear, quark and quantum field theory models. The 
exact energy of the ground state of this oscillator is 
2.044 579 657 447 355 635 36 in atomic units with mass 
set to unity, m = 1. This result is obtained by us by 
a calculation using the Runge-Kutta method in quadru- 
ple precision. The WKB energy is different by 4.5% and 
equals 1.95159 in the same units, while the first-iteration 
QLM energy equals 2.04528 and differs from the exact 
energy only by 0.034%. The QLM energy coincides with 
the exact energy in all twenty digits after the sixth iter- 
ation. 

For the first excited state the exact energy is 
6.713 546 501 445 253 110 53, while the WKB and first- 
iteration QLM energies are 6.656623 and 6.713952 and 
are different from the exact energy by 0.84% and 
0.006% respectively The similar picture exists for 
the second excited state where the exact energy is 
12.767 866 541 180 535 228 88. The WKB and first- 
iteration QLM energies are 12.72396 and 12.76796 and 
arc different from the exact energy by 0.34% and 0.0007% 
respectively. Again, for both first and second excited 
states the QLM energies differ from the exact energies 
only in the twentieth digit after the sixth iteration. 

The graphs corresponding to the Langer WKB solu- 
tion, the exact solution and the first QLM iterate for the 
ground state are displayed in Fig. ^ One can see that 
while the Langer solution is noticeably different from the 
exact solution, the curve of the first QLM iteration is 
indistinguishable from the exact curve. 

This could be followed more precisely by looking at 
Fig. [3 where the logarithm of the difference between the 
exact and WKB solutions and between the exact solution 
and the first QLM iteration are shown. One can see that 
the difference between the exact solution and the first 
QLM iteration is two orders of magnitude smaller than 
the difference between the exact and the WKB solutions, 
that is one QLM iteration increases the accuracy of the 



FIG. 2: Logarithm of the difference between the exact x exact 
and WKB solutions xo (dashed curve) and between the exact 
solution and the first QLM iterate Xm u (solid curve) for the 
ground state of the anharmonic oscillator. The dips on the 
graphs are artifacts of the logarithmic scale, since the loga- 
rithm of the absolute value of the difference of two solutions 
goes to minus infinity at points where the difference changes 
sign. The overall accuracy of the solution can be inferred only 
at x values not too close to the dips. 



result by two orders of magnitude. Note that the dips 
on the graphs are artifacts of the logarithmic scale, since 
the logarithm of the absolute value of the difference of 
two solutions goes to minus infinity at points where the 
difference changes sign. The overall accuracy of the solu- 
tion can be inferred only at x values not too close to the 
dips. 

The accuracy of the WKB approximation increases for 
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FIG. 3: As in Fig.|U but for the excited state of the oscillator 
potential. 



higher excitations. Therefore in the case of the excited 
state both the Langer WKB and QLM curves are indis- 
tinguishable from the exact one. Fig. [3] show, however, 
that also in this case the difference between the exact 
solution and the first QLM iteration is by two orders of 
magnitude smaller than the difference between the exact 
and the WKB solutions. 

Another interesting example is the modified Coulomb 
potential 



V(r) = -- 



1 

■2f> 



h 2 



p 2 (p 



aEr 



which is obtained when the equal masses two-body Dirac 
equation with the static Coulomb interaction is reduced 
to the Schrodinger equation [^H^. The exact energy of 
the ground state with quantum numbers (N, L,S,J) = 
(1,0, 0, 0) is 0.999 993 340 148 538 880 1 in atomic units 
with double mass set to unity, 2 M = 1. This result was 
obtained in the work |22| by an elaborate computation us- 
ing the finite element method and confirmed by ourselves 
using the Rungc-Kutta method in quadruple precision. 
The WKB energy equals 0.999 986 680 and differs from 
the exact one by 6.6 x 10 -4 . The first-iteration QLM en- 
ergy equals 0.999 993 335 and differs from the exact one 
only by 5 x 10 -7 . The QLM energy coincides with the 
exact one in all given digits after the sixth iteration. 

The graph Fig.0]of the exact, WKB and QLM ground 
state wave functions is similar to Fig. [2] 

The graph Fig. [S] for their differences for this case is 
similar to Fig. [21 and shows that the difference between 
the exact wave function and the first QLM iteration is 
by two orders of magnitude smaller than the difference 
between the exact and the WKB solutions. Thus also 
in this case just one QLM iteration increases the accu- 
racy of the wave function by a remarkable two orders of 



magnitude. 

The results for the ground and excited states with 
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FIG. 4: As in Fig.0 but for the ground state with quantum 
numbers (N, L, S, J) = (0, 0, 0, 0) in the modified Coulomb 



potential V(r) = -J- + (1(1 + 1) 
aEr. 
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FIG. 5: As in Fig. [5] but for the ground state in the modified 
Coulomb potential. 



different quantum numbers (TV, L, S, J) for the modified 
Coulomb potential are summed up in Table 1 and also in 
Figs. El |H1 where the the differences between the exact 
wave function and the first QLM iteration and between 
the exact and the WKB solutions are displayed. We see, 
that though the accuracy of the WKB approximation in- 
creases for excited states and states with higher orbital 
momenta, also in these cases one QLM iteration increases 
the accuracy of the wave function by at least two orders 
of magnitude. 
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TABLE I: The WKB, first iteration QLM, full QLM and exact binding energies for different potentials. K denotes the 
number of QLM iterations, m denotes the (reduced) mass of the particle; D\ = 10 2 (£ , oxac t — -E\vKB)/-Eexact, and D2 = 
10 2 (-E CX act — Eq^m) / Eexact ■ (In the first entry, Di and D2 are given in the form x[y] which stands for x x 10 H .) -Ewkb and 

EqI m are given to a limited number of digits, since EqI m is slightly dependent on the discontinuity in the derivative of the 
Langer WKB solution at its joining point between the turning points. This dependence disappears in higher QLM iterations. 
For the high-precision results, the number of correct digits is at least 18, and E tends to be slightly more precise. Consequently 
we give the numbers to 20 or 21 digits. Since the computer arithmetic was quadruple precision (128-bit, about 30 decimal 
places), the differences in the last digits of -Eqlm and -E exac t reflect the different methods used. The state is labeled by nl except 
for the Breit- Coulomb problem where the labels are N, L, S, J. The plus sign for the two-power potential stands for the ground 
(symmetric) state in the corresponding one-dimensional double-well potential; the minus sign stands for the regular states of 
the two-power potential or the antisymmetric state of the corresponding one- dimensional double-well potential, respectively. 
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0.999996670008 
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1101 


0.999996670037 
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0.99999833501727839123 


44 


0.99999833501727839122 


2[-4] 


-8[-8] 






2101 
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1.84744258030447816385 
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Is 
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2.04457965744735563536 
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6.71354650144525311020 
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6.71354650144525311053 


0.85 


-0.006 
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12.72396 


12.76796 


12.7678665411805352297 
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12.7678665411805352289 


0.34 


-0.001 


-24 
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Is 


-17.61192 


-17.5432 


-17.5597967410317970585 
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-17.5597967410317970589 
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0.095 
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-7.190505 


-7.37920 


-7.37854164337449079226 
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-7.37854164337449079262 


2.5 


-0.009 
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-0.029269 
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-0.10819568493119384933 
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-0.009 
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0.49734197 
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0.48314820684089227025 
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0.48314820684089227025 


-2.9 


-0.22 
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1.373747 
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6 


1.37363583606219407958 


-1.5 


-0.008 
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2.17217337 


2.178319 


2. 17745782251542955262 
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2.17745782251542955243 


0.24 


-0.040 



"Includes the tunneling correction to E. 

b Initial WKB approximation includes tunneling correction to E. 



The other examples considered in this paper are the 
logarithmic V(r) = logr, Wood-Saxon V = —Vq/(1 + 
exp((i — R)/a)) and the two-power (double-well) V(r) = 
\g 2 {r 2 — a 2 ) 2 potentials, the results for which are sum- 
marized in Table 1. The graphs corresponding to differ- 
ent states of these potentials are shown in Figs. 151- H7I 
The first two potentials are used respectively for compu- 
tations in quark and nuclear physics. The double-well 
potential, that is the quartic potential in one dimension 
with degenerate minima, is typically studied in quantum 
field theory and in the framework of the tunneling prob- 
lem in quantum mechanics. Its perturbation series does 
not converge and different alternative nonperturbative 
approaches are therefore explored since the description 



of tunneling between two minima should be necessarily 
nonperturbative (see, for example, reference |23j| and the 
references therein). 

In particular, in the paper |24j| using the 1/N expan- 
sion method, the tunneling terms were not included for 
the symmetric (ground) state of the double well poten- 
tial in one dimension, giving the 1/N energy of 0.48305 
compared to the exact energy, 0.48295. . . In addition, in 
our calculation it is easy to specify the boundary condi- 
tion at r = in this particular case (where x(0) 7^ 0), 
so we can calculate on the interval r > only: be- 
cause we do the QLM iteration on the function u(nr) = 
arctan(— Kx( r )/x'( r )); we have simply it(0) = — -|. This 
can easily be seen by taking into account that x( r ) has 
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FIG. 6: As in Fig. |5] but for the excited state with quantum 
numbers (N, L, S, J) = (2, 0, 0, 0) in the modified Coulomb 
potential. 



FIG. 8: As in Fig. |21 but for for the excited state with 
quantum numbers (N, L, S, J) = (2,1,0,1) in the modified 
Coulomb potential. 
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FIG. 7: As in Fig. [5] but for the excited state with quantum 
numbers (N,L,S,J) = (1,1,0,1) in the modified Coulomb 
potential. 



FIG. 9: As in Fig.Q but for the ground state of the logarith- 
mic potential V = log(r). 



an even-power Taylor expansion at r = 0. We use the 
tunneling term just to correct the energy of the initial 
WKB approximation, changing the usual WKB quanti- 
zation condition to 

/ k(r)dr = (n H — )ir e ~5~ a K ( r ) dr 

where the second term on the RHS is the tunneling term 
[25|]; k(r) = iK(r) and n = 0, 1, 2, . . . is the number of 
the bound state. The tunneling correction affects the 1st 
QLM iteration but of course not the full QLM calcula- 
tion, where the boundary conditions completely specify 
the converged solution. 



IV. CONCLUSION 

One can show [TEl Il7j that the approximation by 
the first QLM iterate in Eq. Q leads to exact ener- 
gies for many well known physical potentials such as the 
Coulomb, harmonic oscillator, Poschl- Teller, Hulthcn, 
Hyllcraas, Morse, Eckart, etc. For other potentials which 
have more complicated analytical structure we show on 
examples of the anharmonic oscillator, logarithmic, two- 
power (double- well) , and Wood-Saxon potentials and for 
the solution of the two-body Dirac equation with static 
Coulomb potential, that the use of the Langer WKB wave 
function as an initial guess already in the first QLM ap- 
proximation gives energies and wave functions two orders 
of magnitude more accurate than the WKB results. Such 
a QLM solution, unlike the usual WKB solution, displays 
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FIG. 10: As in Fig. |21 but for the ground state of the loga- 
rithmic potential. 
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FIG. 11: As in Fig. Q but for the ground state of the Wood- 
Saxon potential V = -V /(l + exp((r - R)/a)), with Vo = 24, 
R= 1, a = 0.2. 



no unphysical turning point singularities. Since the first 
QLM iterate is given by an analytic expression Q for 
p = 1 it allows one to estimate analytically the role of 
different parameters and their influence on properties of 
a quantum system with much higher precision than pro- 
vided by the WKB approximation. In addition, it was 
shown that six QLM iterations are typically enough to 
obtain both the wave function and energy with the accu- 
racy of twenty significant digits. 
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FIG. 12: As in Fig. H but for the ground state of the Wood- 
Saxon potential. 
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FIG. 13: As in Fig. but for the first excited state of the 
Wood-Saxon potential. 
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